# Table B.6
rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata0      = paste0(WD,"/_individual_data/_spfrawdata/","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)
library(tidyverse)
library(dotwhisker)
library(DescTools)

load(paste(indata0,"us_spf_pgdp_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data            = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))
data            = ddply(data,.(qdate),transform, fc_rev_avr = mean(fc_rev,na.rm = TRUE))

data$fc_err_pct     = data$fc_err/data$p_realz
data$fc_err_pct_win = Winsorize(data$fc_err_pct, probs = c(0.01, 0.99), na.rm = TRUE)
data$fc_err_win     = data$fc_err_pct_win*data$p_realz

#data$fc_err_win = Winsorize(data$fc_err, probs = c(0.01, 0.99), na.rm = TRUE)

beta_win         =  plm(fc_err_win ~ fc_rev, data=data,index=c("ID", "qdate"), model="within")
double_beta_win  = sqrt(diag(vcovDC(beta_win, type = "HC1")))

delta_win        =  plm(fc_err_win ~ cons, data=data,index=c("ID", "qdate"), model="within")
double_delta_win = sqrt(diag(vcovDC(delta_win, type = "HC1")))

b_win            =  plm(fc_err_win ~ fc_rev_avr, data=data,index=c("ID", "qdate"), model="within")
double_b_win     = sqrt(diag(vcovDC(b_win, type = "HC1")))

stderror = list(double_beta_win, double_b_win, double_delta_win)

stargazer(beta_win, b_win, delta_win, 
          type = "text",
          se = stderror,
          out="_tables/table_b6.txt")